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ABSTRACT 

This chapter describes recent developments which have improved our under- 
standing of how finite difference methods resolve discontinuous solutions to hyper- 
bolic partial differential equations. As a result of this understanding improved shock 
capturing methods are currently being developed and tested. Some of these meth- 
ods are described and numerical results are presented showing their performance 
on problems containing shocks in one and two dimensions. 

We begin this discussion by defining what is meant by a conservative differ- 
ence scheme and showing that conservation implies that, except in very special 
circumstances, shocks must be spread over at least two grid intervals. These two 
interval shocks are actually attained in one dimension if the shock is steady and an 
upwind scheme is used. By analyzing this case, we determine the reason for this 
excellent shock resolution and use this result to provide a mechanism for improving 
the resolution of two dimensional steady shocks. Unfortunately, this same analysis 
shows that these results cannot be extended to shocks which move relative to the 
computing grid. 

To deal with moving shocks and contact discontinuities we introduce total 
variation diminishing (TVD) finite difference schemes and flux limiters. We show 
that TVD schemes are not necessarily upwind, but that upwind TVD schemes 
perform better because they permit a wider choice of flux limiters. The advantage 
of non-upwind TVD schemes is that they are easy to implement. Indeed, it is 
possible to add an appropriately chosen artificial viscosity to a conventional scheme 
such as MacCormack’s method and make it TVD. We conclude by presenting some 
theoretical results on flux limiters and some numerical computations to illustrate 
the theory. 


Research was supported by the National Aeronautics and Space Administration 
under NASA Contract No. NASl-17070 while the author was in residence at ICASE, 
NASA Langley Research Center, Hampton, VA 23665. 





1. Introduction 


We wish to consider the numerical solution of hyperbolic systems of the form 

U( + /(u)a; =0, f > 0, -00 < X < 00 (1.1) 

where u(a:,<) and /(u(x, t)) are m-vectors. 

In general, systems of this form admit discontinuous (shock) solutions which 
satisfy the Rankin e-Hu goniot jump conditions 

S|u 1=1/1 (1.2) 


where | u | denotes the jump in the variable u across the shock and s is the speed 
of propagation of the shock. 

In the following we consider finite difference schemes which can be put into the 
following conservation form 

--^”-1/2 W-i-^ (1-3) 

because Lax and Wendroff I"*! have shown that covergent schemes of the form (1.3) 
converge to solutions which satisfy equation (1.2). To show why these schemes are 
called conservative, we apply (5.2.3) to (1.1) with periodic boundary conditions and 
sum over one period. The result is 

= ; (1.4) 


if we apply this result recursively, we obtain 

E = E 


which says that the discrete integral of the solution is conserved. This conserva- 
tion property implies that shocks which fall between mesh points can at best be 
represented by transitions over two intervals. That is, a shock of the form 


u(x) = I 


for X < Xj 
for X > X* 


where JCj- 1/2 < x® < Xj^i /2 can, at best, be represented by a transition of the 
form 

{ U£, for I < j 

Urn, iovi = j (1.7) 

Ufl, for t > j 
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where 


“m (xy+i/2 - Xy_i/2) = Ul(x* - ~ U/i(xy+i/2 - X*) 

Xj+1/2 = (Xj+1 +Xj)/2 

satisfies (1.5). In most practical computations, shocks are spread over many more 
than two intervals. In this paper we examine why this is so and how shock resolution 
approaching that of equation (1.7) can be attained. 

2. One Dimensional Steady Shocks 


In this section we study the application of upwind diflFerence schemes to the 
numerical computation of one dimensional steady shocks. We undertake this study 
not because this problem is of general interest but because it is known that upwind 
differences schemes attain two interval shock resolution in this case. By studying 
this case in detail we can determine which features of these schemes are responsible 
for this excellent shock capturing ability. In the next section we use this knowledge 
to construct shock capturing schemes for two dimensional steady shocks. Since this 
theory requires moving grids to capture moving shocks, we examine the moving 
shock problem from a difiFerent point of view in a later section. 

To simplify the details we restrict this discussion to the scalar equation 


«f + f{u)x = ut + a(u)«a = 0, a(u) 


— t > 0, Xq < X < Xjv 
du 


( 2 . 1 ) 


and consider the first order upwind difference scheme of Murman and Cole This 
is a reasonable thing to do because modern high order methods are designed to 
become first order in the vicinity of shocks and because the Murman and Cole 
method can be extended to systems by a procedure developed by Roe 


Consider a partition of the x axis 



Xo < a:i < ... < Xat 

(2.2) 

and locally linearize equation (2.1) on each interval (xy, xy+i). That is let 



Ut +a(uy, Uj+i)Ux = 0 

(2.3) 

where 

a(u ■ u 0 . 1 1 = 1 ^ ® 

1 j+iJ ydf[uj)/du, if A«y = 0 

(2.4) 

and A«j = «,+i 

— «f. A/,- = fi+i — ft. If we define 



a‘^(«y, Uj+i) = max (O,o(uy, Uy+i)) 

(2.5) 

and 

a~(uj, uy+i) = min (0, o(«y, «y+i)) , 

(2.6) 
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then we can write the Murman-Cole scheme in the form 


jjn+l 


Uf 


2At 




a+{U^_^,U^)AU;f., 


( 2 . 7 ) 


+ a-{Uf,U‘+,)AUf 


or in the alternative form 




This method can also be written in conservation form 

2At 


(A»,„ + A»,) {nUhUf^Cl-F(V;_,,Un)} 


( 2 . 8 ) 


( 2 . 9 ) 


if we define 


and note that 


F(Uj, «y+i) = /(<lj+l) - A/t 
= /(«,) + A/- 

= l/2{/(ay) + + AfJ - A/t} 


/(")+i) - /(%) = A/t + A/j. . 


( 2 . 10 ) 


( 2 . 11 ) 


The main result of this section is the following theorem. 

Theorem . Consider the following representation of a steady shock at time V 



for y < 0 


^771) 

for j = 0 

(2. 


for y > 0 


where Um is some intermediate value between and ur while and ur satisfy 
the steady Rankiue-Hugouiot jump condition 


f{uR) - /(ui) = 0 

and the geometric entropy condition 




if Um is such that 


o(«t,«m) >0>a(Um,!lfi) 
theu (2.12) is a stable steady solution to (2.S). 


( 2 . 13 ) 

( 2 . 14 ) 

( 2 . 15 ) 
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Proof A direct calculation shows that equation (2.12) solves (2.8). That b: 
for j < —1 

A/t = /(ui) - /(ul) = 0 

AfJ' = 0 , since a{uL, ul) = o(u£,) > 0 

for y > 1 

A// = 0 , since a(ufl, ur) = o(«/j) < 0 
A/- = /{ur) - f{uR) = 0 

for j = —1 

A/^i = f{Um) ~ /(wl) 

A/ri = 0 , since a{uL, Um) > 0 

and for j = 0 

A/^ = 0 , since a(«m) «r) < 0 
A/o~ = /{ur) - f{um) 
so 

A/j?! + A/q = f{Um) - /(«l) + /(«k) - /(«m) 

= /{ur) - /(U£,) = 0. 

To prove stability assume that u'^ = ui, + e''^ for some j < —1. Then 


{ul + e”+^) = {uL+e^)- -^[/(«i + e*") - /(«£,)] 
= Ui + e“ - ^la(uL)e" + 0(e”)“| 

= (l-«L|^)e”+0(e”f. 


Since 

At , At^ 

gw+i 0 as n — > 00 and the solution is stable. 

A similar argument holds for i > 1. | 

If we consider the conservation form of the scheme (2.9), (2.10), we note that 
the numerical fluxes 

F{UL,Um) = f{UL) (2.16) 

and 

F{Um,UR) = f{uR) (2.17) 

do not depend on the intermediate value Um- If we study the pointwise error in 
the computed solution, we note that since Um should have the value of either 
or UR,Um is in error by order («r — «£,). Equations (2.16), (2.17) trap this error 
at one grid point and prevent it from propagating through the domain. This is the 
property of the upwind methods which is responsible for their excellent ability to 
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capture steady shocks. Methods with a numerical flux that depend explicitly on Um 
will create an error that will propagate throughout the domain of the calculation. 
This error will show up either as a spreading of the shock or as wiggles in the 
solution. 

Unfortunately, if a shock moves through a grid, the computed solution at points 
left behind will depend on intermediate values such as Thus, it is not possible 
to avoid the propagation of error to points away from a shock when a shock moves 
relative to the computing grid. 

3. Two Dimensional Steady Shocks 

In the preceding section we showed how upwind differencing improves the res- 
olution of steady one dimensional shocks. Here we extend these results to the two 
dimensional Euler equations of gas dynamics. 

To show that more is needed than a straightforward application of the results 
of the last section, we consider the exact weak solution to the Euler equations 
which corresponds to a plane oblique shock. Such a solution can be constructed 
by superimposing an arbitrary tangential velocity on a one dimensional shock and 
choosing a cartesian coordinate system that is aligned with the upstream velocity 
vector (see Figure l). 



Figure 1. Plane oblique shock solution. 

Before we continue, we note that for the Euler equations the direction of differ- 
encing for the X and y derivatives in an upwind scheme is determined by the local 
Mach numbers and My, respectively, where 


Mx = u/a, My = v/a, a = y/'^pfp 


(3.1) 




and («, w) are the components of the velocity vector v along the components (a:, y) of 
the position vector x. If Mx > 1, backward differences are used for all x derivatives. 
If Mx < —1, forward differences are used for all x derivatives. If —1 < Mx < l,x 
derivatives are approximated by both forward and backward differences according 
to some recipe (cf., e.g., Osher^^^, Roe^®], van Leer^^^J). Corresponding rules hold 
for the y derivatives. 

Returning to the oblique shock solution illustrated in Figure 1, we note that 
since the tangential velocity is arbitrary, it can be chosen large enough to make 
Mx > 1 on both sides of the shock. If this is done, the direction of differencing does 
not change across the shock, as it should. The method does not “see” the shock 
and the shock is spread as shown in Figure 2. 

Other difficulties involve the differencing of derivatives in the direction tangent 
to the shock. We discuss these problems below. 

One way to avoid these problems is to choose a computing grid which is aligned 
with the shock. This restores the one dimensional shock resolution. Unfortunately, 
it is very difficult to choose such a grid, especially when multiple shocks intersect. 
For this reason we pursue the alternative approach of deriving a method on an 
arbitrary grid which mimics the behavior of a method constructed on an aligned 
grid. 

To carry this out, we write the Euler equations in a local coordinate system, 
{x',y') which is aligned with potential shock directions and is rotated with respect 
to the coordinate system of the computing grid (*, y). At this time we assume that 
these directions are known. This geometry is shown in Figure 3. 


Since the Euler equations are invariant under rotation, they can be written 
immediately in the new coordinate system as 


dW _ df'{W) _ dg'{W) 

dt dx' dy' 

where 

w' = [p,pu',pv',e\^ 
f' = [pu', pu'^ + p, pu'v', (e + p)«']^ 
g' = \pv\ pu'v' , pv'^ + p, (e + p)v']^ 
e = p/(7 “ 1) + l/2p(u'^ + 

and 

u' = u cos 0 + V sin 0 
v' = — u sin 0 + V cos 0. 


(3.2) 


(3.3) 


(3.4) 


If the coordinate system is aligned with a plane steady shock, the derivatives in 
the second term on the right of (3.2) would exist and be equal to zero. This reduces 
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the two-dimensional problem, locally, to the one-dimensional problem described in 
the previous section. By choosing different discretizations for each term on the right 
of (3.2), we attempt to construct a numerical scheme which mimics this behavior. 
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Figure 2. Density profile through an oblique shock computed using a first 
order upwind scheme. 



Figure 3. Geometry of local and global coordinate systems. 
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If we express the first term on the right of (3.2) in global coordinates, we get 

(3.6) 


3f' . dt' 

_ = cos«^ + sm»— . 


This is approximated by the finite difference expression 

- F'(WUi^Wr)\ 

+"^[F'(Wl,, (3.6) 

where F' is an upwind numerical fiux such as that described in section 2.3 and W 
is the numerical approximation to w'. 


The second term on the right of (3.2) is the directional derivative of g' in the 
y' direction. If we write this expression in global coordinates we get 


— = -sin0— +COS0— — 
ay' dx dy 


(3.7) 


when we attempted to approximate (3.7) on the usual five point stencil, the results 
were unacceptable regardless of the choice of numerical fiux. This leads us to 
believe that an acceptable approximation to (3.7) must include the corner points 
(®t+i)%+i) 3^d (ajj-i, j/y_i). Therefore, we chose the discretization 


+ 


Aa: 

COS0 


Ay 


where the functions iT and *2 are chosen so that the resulting stencil provides 

the closest possible approximation to the directional derivative. These functions are 
tabulated in Table 5.1 and the resulting stencils are shown in Figure 4. With this 
construction, the results did not seem to depend strongly on the choice of numerical 
flux although less dissipative numerical fluxes gave slightly better results than more 
dissipative numerical fluxes. 


Case 

Angle Range 

jl 

j2 

il 

i2 

a 

-00 < tan(0) < -AxjAy 

j 

j 

i-1 

i+1 

b 

—AxjAy < tan{9) < 0 

M 

j+1 

i 

i 

c 

0 < tan{6) < AxjAy 

j+i 

M 

i 

i 

d 

AxjAy < tan{6) < oo 

j 

j 

i+1 

i-1 


Table 5.1 Definition of Computational Stencils. 
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Figure 4 Computational stencils for tangential derivative approximations. 

Expressions (3.6) and (3.8) have been derived under the assumption that we can 
associate a unique value of the angle Q with each computational cell. Unfortunately, 
it does not seem to be possible to allocate the angles in this way and still construct 
a conservative difference scheme. Therefore we have chosen to associate a value 
of the angle with each cell boundary. This permits us to construct the following 
conservative difference scheme 


nr = n.i - s - fi [G,j+,/, - G( 


(3.9) 


where cartesian tensor transformation rules are used to construct the numerical 
fluxes from the rotated numerical fluxes. That is, formula of the form 


fi[w) = pu = pu' cos 6 — pv' sin 6 

— f'i{w')cosO — gi(u;')sin 0 


(3.10) 


is used when fi is the flux of a scalar variable such as p or e and a formula of the 
form 

/2 (ty) = +p = p[u' cos 9 — v' sin 0 )^ + p 

= pu'^ cos^ 9 — 2pu'v' cos 0 sin 0 + pv'^ sin^ 9 + p(cos^ 9 + sin^ 9) (3-H) 

= ^ ~ cos 0 sin 9 — { 73 ( 1 ^') cos 0 sin 0 + /j(w')sin^ 9 

is used when /2 is the flux of a vector component such as pu or pv. 
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Finally, we show how the angles used in the algorithm are determined. We note 
in passing that we do not attempt to determine whether or not a shock actually 
exists. Instead the algorithm assumes that a shock exists between each pair of grid 
points and determines its orientation. This approach locates the proper direction 
when it is needed and causes no problems otherwise. 

In reference [l], we determine the direction of the normal to a shock which 
is assumed to lie between two grid points as the direction of the jump in velocity 
between the grid points. This works well when the computing grid is uniform but 
fails when the grid is highly stretched. To overcome this problem we now determine 
the direction of a shock normal as the direction of the gradient of a scalar variable 
such as the pressure or density. Suitably smoothed numerical approximations to 
the gradient have provided proper shock angles in all cases that we have tried. 

Figures 5 to 7 show the results of computations of a two dimensional shock 
reflection using the first order scheme of Osher, the second order scheme of van 
Leer, and the present scheme respectively. These figures clearly show that rotational 
bias considerably improves the shock resolution of first order schemes. Indeed the 
shock resolution of a first order accurate rotationally biased scheme rivals that of a 
sophisticated second order method. 

In a first attempt to obtain a higher order rotational scheme, we replaced 
the first order numerical fluxes of the present method with second order numeri- 
cal fluxes. The resulting scheme gave some improvement in shock resolution but 
was quite sensitive to the treatment of boundary conditions and was probably too 
complex to be of practical use. It is not clear that the small improvement in shock 
resolution over current second order methods would justify the considerable increase 
in complexity of the scheme. It is also not entirely clear that the resulting scheme 
is actually second order. For these reasons there is still need for research on high 
order two dimensional schemes. 

4. TVD Finite Difference Schemes and Flux Limiters 

The method described in the previous section is at this time only first order 
accurate and has been designed especially for problems with steady shocks. In this 
section, we examine second order methods which are designed to give good results 
when applied to problems containing moving shocks and contact discontinuities as 
well as steady shocks. 

In order to demonstrate the main ideas in a simple setting, we consider here 
only the linear advection equation 

ut + aux =0, a = constant. (4.1) 

The extension of this work to nonlinear equations, systems and multiple dimensions 
is described in reference [2]. 
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We consider explicit finite difference schemes in conservation form which ap- 
proximate equation (4.1) and which we denote by 

= L.U^. (4.2) 

The total variation of a mesh function is defined by the formula 

TV(V]=Y.Pl^.,-Uf\ (4.3) 

t 
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I 


Figure 5 Three dimensional density 
plot for shock reflection computed 
order rotated difference method of 
Davis 


Figure 6 Three dimensional density 
plot for shock reflection computed 
using the second order method of 
van Leer 


Figure 7 Three dimensional density 
plot for shock reflection computed using 
the first order method of Osher and 
Solomon 



and a scheme is called total variation diminishing (TVD) if 

2V(?7''+i) = TV[L.U^) < TV{U^). (4.4) 

In the following we consider only total variation diminishing schemes because solu- 
tions to TVD schemes do not exhibit spurious oscillations. 

To determine whether a scheme is TVD, we rewrite it in the form 

t'r* = - Ot-,/2 (4.5) 

where 

- n (4.6) 

and C'fc_i/ 2 >-Dfc+i /2 2 tre functions of and we apply the following result. 

Lemma (Harten). If the coefficients G and D of equation (5.4.5) satisfy the in- 
equalities 

0 ^ <^*+ 1/2 

0 < I>fc+i/2 (4.7) 

0 < <^*+ 1/2 + ■^^fc+ 1/2 ^ 

then the scheme (4.5) is TVD. | 

There are a number of ways to construct TVD schemes. See the papers by 

Harten Osher^^J, Roe^^l, and van Leert^^J for details of other approaches. Here 
we construct a scheme by adding a term of the form 


■^*+l/2^^fc+l/2 ~ -^A-l/2^^fc— 1/2 (^•®) 

to the Lax-WendroflF method 

= n - (4.9) 

where u = aAtfAx and K in equation (4.8) is chosen so that the resulting scheme 
is TVD. 

If a > 0 and K is chosen to have the form 


where 


Kk+i /2 = - 4>{rt)] 


+ 

k 




(4.10) 


(4.11) 
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and <j> is the flux limiter, we obtain a scheme of the form (4.5) studied by Swebyl^^l 


C'fc-i/2 = i^{l+ 1/2(1 - i^)[(l>{rt)/rt - 0(r^_i)]} 

Dk+i = 0 . 


This scheme will be TVD provided u < 1 and 

-2 


1 - 1 / 




(4.12) 


(4.13) 


If a < 0 and K is chosen to have the form 

1/(1 + u) 




where 




= 


we obtain a scheme of the form (4.5) with 

^k-il2 — 0 

Dk+i /2 = -1 + 1/2(1 + v) [0(r-+i) - 0(r^ )/r“ ] }. 

This scheme will be TVD provided —u < 1 and 

_ z 2 l . < _ 0(rr. J1 < — . 

I + I' “ ^ rr ^ “ V 


(4.14) 


(4.15) 


(4.16) 


(4.17) 


The schemes (4.12) and (4.16) can be combined into a single upwind scheme if 
K is chosen to have the form 


/ V 


K, 


fc+l/2 


= < 


f(l-z/) [l-«^(rj^)l , ifa>0 

1(1 + 1/) [<?i(r^+i)-l], ifo <0 


(4.18) 


For hyperbolic systems it would be convenient to have a method that did not require 
that we know which direction is upwind. To this end we rewrite (4.18) in the form 

Kk+if2 = y(l - IH)[1 (4.19) 

where <j>' is an appropriately determined flux limiter. Roel^^l has done a considerable 
amount of work on what constitutes an appropriate flux limiter. In the remainder 
of this section we discuss this work. 
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In order to assure that the scheme (4.8), (4.9) is TVD, we require that inequal- 
itites (4.13) and (4.17) hold. To simplify the analysis, we also require that <f) > 0 
and 0 = 0 when r < 0. This reduces the scheme to first order in the vicinity of 
extrema and causes peaks to be “clipped.” Attempts to correct these problems are 
the topics of current research. With these simplifications, inequalities (4.13) and 
(4.17) become 


(r^) < min 


_ 2 _' 
u ’ 1 - I/. 


(4.20) 


and 


respectively. 




< 


min 


' 2r^ 2 • 

. IH ’ 1- M. 


(4.21) 


Roe has shown that the scheme will be second order accurate if 


0(r) = r + 0(1 - r) (4.22) 

in the vicinity of r = 1. There is also a condition which assures that the method 
will be third order accurate, but it requires that the upwind direction be known so 
we will not discuss it here. 


Of more interest is the limiter 


0(r) = 


2r ,(l — r*^) — 1/(1 — r) 
u{l - u) (1 - r)2 


(4.23) 


derived by Roe which exactly convects exponentially varying data. Although this 
limiter which Roe calls Hyperbee, is too complex and does not perform well in 
practice, it seems to provide a yardstick for choosing good limiters. In particular. 
Roe’s numerical experiments seem to show that any limiter which crosses Hyperbee 
in two places other than r = 1 will convect step data as a smooth profile which does 
not spread with time. This is an exciting result since all other schemes are known 
to spread step data. 


A closer look at these results shows that this smooth profile consists of two 
exponentials, corresponding to those values of r where the limiter crosses Hyperbee, 
connected by a smooth curve that does not grow in time. Recently, this author has 
constructed mathematical proofs, of some of these results. 


Figure 8 is a sketch of the Hyperbee limiter and the upwind limiter 


0(r) = max 


0, mm 



min 


r, 


2 

1-H 


)i 


(4,24) 


which Roe calls Ultrabee. Figure 9 shows the results of a computation using the 
Ultrabee limiter to convect a square wave for 250 time steps at a Courant number 
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of 0.5. Figure 10 shows the results of this same computation using the non-upwind 
limiter 


2?*”^ 2 2y*~” 2 

0'(r+, r") = max{0, r“, iHH 

These figures show that the upwind limiter gives sharper steps but the non-upwind 
limiter still gives acceptable results. The upwind limiter is able to satisfy more 
closely the inequalitites (4.20), (4.21) and this is the reason for its better perfor- 
mance. On the other hand, the non-upwind limiter is simple to program and would 
probably run faster than the upwind limiter. 



Figure 8. Flux limiters. 




Figure 9. Linear advection of a square pulse computed using the MacCormack 
scheme and the upwind Ultrabee limiter. 



Figure 10. Linear advection of a sqaure pulse computed using MacCormack scheme 
and the non-upwind limiter equation (4.25). 
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X 


Figure 11. Three dimensional density plot for shock reflection computed using the 
MacGormack scheme and the non-upwind limiter equation (4.25). 

Figure 11 shows results of the application of the non-upwind limiter to the 
shock reflection problem. We believe that the standing waves along the reflected 
shock are due to the treatment of the wall boundary conditions or to the way that 
we extended our one dimensional results to two dimensions. This is currently being 
investigated. 


5.6 Summary 

In this chapter we have examined how upwind differencing and flux limiters 
improve the resolution of numerical approximations to discontinuous solutions of 
hyperbolic systems. The study of these numerical processes has led to the devel- 
opment of new numerical methods of high accuracy which resolve shocks without 
excessive spreading or spurious wiggles. 

Here we have concentrated on three special cases. In particular we showed 
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how the requirement that a finite diflference scheme be conservative is equivalent to 
the requirement that some pointwise error be made in the vicinity of a shock. On 
the other hand a properly constructed upwind difference scheme will confine these 
errors to one mesh interval if the shock is steady. This result has been applied to 
one dimensional problems for some time. Here we also showed how to apply it to 
two dimensional problems. 

For moving discontinuties, it is well known that the use of flux limiters permits 
the construction of high order upwind schemes which can resolve these discontinu- 
ities without spurious wiggles. Here we have shown that this construction can also 
be applied to second order central difference schemes and that the resulting schemes 
are simpler than their upwind counterparts. 

Finally, we have described an exciting discovery by P. L. Roe. Roe has charac- 
terized a class of flux limiters which are observed experimentally to approximate a 
discontinuity by a smooth narrow transition that does not spread as it is advected. 
So far as we know all other schemes spread linear discontinuities for all time. 
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